home *** CD-ROM | disk | FTP | other *** search
/ The Datafile PD-CD 5 / DATAFILE_PDCD5.iso / utilities / m / mole / PDB->Mole (.txt) < prev   
RISC OS BBC BASIC V Source  |  1993-04-21  |  6KB  |  214 lines
  1.  PDB to Mole file converter
  2.  Simon Kilvington, 1993
  3.  see "Protein Data Bank Atomic Coordinate Entry Format Description"
  4.  published by Brookhaven National Laboratory, for a description of
  5.  the PDB file format
  6.      version$="v0.02 (21-Apr-93)"
  7. @maxatoms=512:
  8.  you can adjust this if you have enough memory
  9.  atomname$(50),connectdist(50),atmrad(50),atmcol%(50)
  10.  serial%(maxatoms),name$(maxatoms),residue$(maxatoms),resseq%(maxatoms)
  11.  xpos(maxatoms),ypos(maxatoms),zpos(maxatoms)
  12.  connect%(maxatoms),connectto%(maxatoms,3)
  13.  read in connection distance data
  14. Edisterr=1:
  15.  error in coordinates (used later in connecting atoms)
  16. ("PDBFiles.connect")
  17. line$=
  18. getline(H%)
  19.  line$<>"" 
  20. !atomname$(A%)=
  21. getword(line$)
  22. (atomname$(A%))=1 atomname$(A%)=" "+atomname$(A%)
  23. :connectdist(A%)=
  24. getword(line$))*100:
  25.  convert to pm
  26. !atmrad(A%)=
  27. getword(line$))
  28. +atmcol%(A%)=
  29. ("&"+
  30. getword(line$)+"00")
  31.     A%+=1
  32. lasttype%=A%
  33.  get the input and ouput file names
  34. "PDB to Mole file converter "+version$
  35. "==========================="+
  36. (version$),"=")
  37. '"PDB files are taken from directory 'PDBFiles.pdb'"
  38. "Mole files are saved in directory 'MoleFiles'"
  39. '"*** WARNING *** dont use this to convert large protein files as it is not finished yet!"
  40. '"PDB file leafname  : "pdbfile$
  41. "Mole file leafname : "molefile$
  42. '"Reading in PDB data..."
  43. 2"H%=
  44. ("PDBFiles.pdb."+pdbfile$)
  45. A%=0:C%=0
  46. record$=
  47. getline(H%)
  48. record$,6) 
  49.  "ATOM  ","HETATM"
  50. 8 serial%(A%)=
  51. record$,7,5))
  52. name$(A%)=
  53. record$,13,4)
  54. residue$(A%)=
  55. record$,18,3)
  56. ;!resseq%(A%)=
  57. record$,23,4))
  58. <2xpos(A%)=
  59. record$,31,8))*100:
  60.  convert to pm
  61. ="ypos(A%)=
  62. record$,39,8))*100
  63. >"zpos(A%)=
  64. record$,47,8))*100
  65. ?    A%+=1
  66.  "SSBOND"
  67. "Record type 'SSBOND' not yet supported, so ignored"
  68.  "CONECT"
  69. C!connect%(C%)=
  70. record$,7,5))
  71. T%=0 
  72. record$,12+T%*5,5)
  73.  to$="     " connectto%(C%,T%)=-1 
  74.  connectto%(C%,T%)=
  75. (to$)
  76. H    C%+=1
  77.  "TER   "
  78. "Record type 'TER   ' not yet supported, so ignored"
  79.  "REMARK","FTNOTE"
  80.  record$
  81. "Unknown record type '";
  82. record$,6);"' ignored"
  83. lastatom%=A%
  84. lastconnect%=C%
  85. '"Building Mole file..."
  86. decpt=8
  87. hdrsize%=4+4+4+8+4+12+12+4
  88. atmsize%=12+4+4+4+16
  89.  data% hdrsize%+lastatom%*atmsize%
  90. dataend%=data%
  91.  set up a default header
  92. ^>!dataend%=&454C4F4D:dataend%!4=200:
  93.  id word, file version
  94. _,dataend%!8=&00424450:
  95.  title 'PDB'+CHR$0
  96. `Bdataend%!12=&00000000:dataend%!16=&00BBFF00:
  97.  bg and bond cols
  98. a8dataend%!20=1<<decpt:
  99.  scaling factor (1 OS unit/pm)
  100. b8dataend%!24=0:dataend%!28=0:dataend%!32=0:
  101.  rotation
  102. c;dataend%!36=0:dataend%!40=0:dataend%!44=0:
  103.  translation
  104. d#dataend%!48=lastatom%:
  105.  # atoms
  106. dataend%+=hdrsize%
  107.  now add the atom data
  108.  "atom ";A%+1;" of ";lastatom%
  109. i" !dataend%=xpos(A%)*(1<<decpt)
  110. j# dataend%!4=ypos(A%)*(1<<decpt)
  111. k# dataend%!8=zpos(A%)*(1<<decpt)
  112.  get the part of the atom name to check against our atom types list
  113.  symbol1$=
  114. name$(A%),1)
  115.  symbol2$=
  116. name$(A%),2,1)
  117.  symbol2$>="A" 
  118.  symbol2$<="Z" symbol$=" "+symbol2$ 
  119.  symbol$=symbol1$+symbol2$
  120. p     B%=0
  121.  atomname$(B%)<>symbol$
  122.   B%+=1
  123.  B%=lasttype% 
  124. "Unknown atom type '"+symbol$+"'":
  125. u& dataend%!12=atmrad(B%)*(1<<decpt)
  126.  dataend%!16=atmcol%(B%)
  127.  calculate the valency from what we know
  128. x     C%=0
  129.  C%<lastconnect% 
  130.  connect%(C%)<>serial%(A%)
  131.   C%+=1
  132.  val%=0
  133.  C%<lastconnect% 
  134.  connection data was in the file
  135. V%=0 
  136.  connectto%(C%,V%)<>-1 val%+=1
  137.  see if we have a template for this residue type
  138. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  139.  this inability to use connectivity templates means it will take an incredibly long time to
  140.  convert large protein data files, and even when it does eventually finish it may well have
  141.  added extra bonds to the model
  142. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  143.  no template so use the distance data to see whats connected to what
  144. atm%=0 
  145.  lastatom%-1
  146.  atm%<>A% 
  147.      symbol1$=
  148. name$(atm%),1)
  149. "    symbol2$=
  150. name$(atm%),2,1)
  151. X    
  152.  symbol2$>="A" 
  153.  symbol2$<="Z" atmsym$=" "+symbol2$ 
  154.  atmsym$=symbol1$+symbol2$
  155.     R%=0
  156. !    
  157.  atomname$(R%)<>atmsym$
  158.      R%+=1
  159. <     
  160.  R%=lasttype% 
  161. "Unknown atom type '"+atmsym$+"'":
  162.         
  163.     dx=xpos(atm%)-xpos(A%)
  164.     dy=ypos(atm%)-ypos(A%)
  165.     dz=zpos(atm%)-zpos(A%)
  166. N    
  167. (dx*dx+dy*dy+dz*dz) < (connectdist(B%)+connectdist(R%)+2*disterr) 
  168. *     connectto%(C%,val%)=serial%(atm%)
  169.      val%+=1
  170.         
  171.  val%=0 
  172.  val%>4 
  173. "Atom ";A%;" is unconnected or over connected":
  174. 0 dataend%!20=((val%-1)<<16)+val%:
  175.  info word
  176.  symbol$=" H" dataend%!20+=&80000000
  177.  dataend%+=24
  178.      B%=0
  179. V%=1 
  180.  val%
  181.  connectto%(C%,B%)=-1
  182.    B%+=1
  183.   atom%=0
  184.  serial%(atom%)<>connectto%(C%,B%)
  185.    atom%+=1
  186.  atom%=lastatom% 
  187. "Atom ";A%;" is connected to a non-existant atom":
  188. !  !dataend%=atom%:dataend%+=4
  189.   B%+=1
  190.  A%+=1
  191.  A%=lastatom%
  192. '"Saving as 'MoleFiles."+molefile$+"'"
  193. ("Save MoleFiles."+molefile$+" "+
  194. ~data%+" "+
  195. ~dataend%)
  196. ("SetType MoleFiles."+molefile$+" Mole")
  197. '"You may need to recentre the molecule before you can see it in Mole"
  198. getline(H%)
  199. line$=""
  200. byte%=
  201.  byte%>=32
  202. line$+=
  203. byte%
  204. byte%=
  205. =line$
  206. getword(
  207. T$,1)=" "
  208. T$,2)
  209.     W$=""
  210. (T$)>0 
  211. T$,1)<>" "
  212. T$,1)
  213. T$,2)
  214.